#!/bin/bash
set -e

function info() {
echo Usage: `basename $0` [-q mapping_quality] in.bam in.bed
exit 65
}

while getopts  ":q:p:" opts
do
        case  $opts  in
        q) quality=$OPTARG;;
		p) out_prefix=$OPTARG;;
		*) info;;
        esac
done
shift $(($OPTIND - 1))


if [ $# -lt 2 ]; then info; fi


. /mnt/ilustre/app/medical/tools/.var

if test -z $quality; then quality=1; fi

interval0=$out_prefix.bed
# cut -f1-3 $2 > $interval0
cat $2 >$interval0

echo;echo samtools view mapping quality '>=' 1 and no dup
samtools view -uF 0x400 -q$quality $1| coverageBed -abam - -b $interval0 -d >& $out_prefix.cov0.txt

echo;echo cal cov
cal_coverage_metrics.pl $out_prefix.cov0.txt $out_prefix.cov.txt

